Methods and apparatuses for estimating the elliptical cone of uncertainty

ABSTRACT

A magnetic resonance imaging scanner acquires diffusion-weighted imaging data. A reconstruction engine reconstructs the acquired diffusion-weighted imaging data into diffusion-weighted image representations. A diffusion tensor engine constructs a diffusion tensor map of an area of an interest of a subject. An eigenvalue/eigenvector ordering engine obtains and orders eigenvectors and eigenvalues at each voxel. A covariance matrix determining engine constructs a covariance matrix of a major eigenvector of each voxel. A first normalized measure determining engine computes a first normalized measure. A second normalized measure determining engine computes a second normalized measure. A rendering engine generates a human-viewable display of an image representation.

This application claims the benefit of U.S. Provisional Application No.60/996,169 filed on Nov. 5, 2007, of the common assignee, entitled“Methods and Apparatuses for Estimating the Elliptical Cone ofUncertainty,” the contents of which are incorporated by reference intheir entirety.

BACKGROUND

The following relates to the magnetic resonance arts. It particularlyrelates to the imaging, tracking, and displaying of neural fibers andfiber bundles by diffusion tensor magnetic resonance imaging (DT-MRI),and will be described with particular reference thereto. However, itwill also find application in conjunction with the tracking andgraphical rendering of other types of fibrous structures as well as withother imaging modalities such as single photon emission computedtomography imaging (SPECT), computed tomography (CT), positron emissiontomography (PET), and also in crystallography and geophysics.

Nerve tissue in human beings and other mammals includes neurons withelongated axonal portions arranged to form neural fibers or fiberbundles, along which electrochemical signals are transmitted. In thebrain, for example, functional areas defined by very high neuraldensities are typically linked by structurally complex neural networksof axonal fiber bundles. The axonal fiber bundles and other fibrousmaterial are substantially surrounded by other tissue.

Diagnosis of neural diseases, planning for brain surgery, and otherneurologically related clinical activities as well as research studieson brain functioning may benefit from non-invasive imaging and trackingof the axonal fibers and fiber bundles. In particular, diffusion tensormagnetic resonance imaging (DT-MRI) has been shown to provide imagecontrast that correlates with axonal fiber bundles. In the DT-MRItechnique, diffusion-sensitizing magnetic field gradients are applied inthe excitation/imaging sequence so that the magnetic resonance imagesinclude contrast related to the diffusion of water or other fluidmolecules. By applying the diffusion gradients in selected directionsduring the excitation/imaging sequence, diffusion weighted images areacquired from which diffusion tensor coefficients are obtained for eachvoxel location in image space.

The intensity of individual voxels are fitted to calculate sixindependent variables in a 3×3 diffusion tensor. The diffusion tensor isthen diagonalized to obtain three eigenvalues and corresponding threeeigenvectors. The eigenvectors represent the local direction of thebrain fiber structure at the voxel at issue. Within an imaging voxel,the directional information of white matter tracts is usually obtainedfrom the major eigenvector of the diffusion tensor.

Using the local directional information, a global anatomicalconnectivity of white matter tracts is constructed. Each time thetracking connects a pixel to the next pixel, judgment is made whetherthe fiber is continuous or terminated based on randomness of the fiberorientation of the adjacent pixels. The interplay between the localdirectional and global structural information is crucial inunderstanding changes in white matter tracts. However, image noise mayproduce errors in the calculated tensor, and, hence, in its eigenvaluesand eigenvectors. There is an uncertainty associated with every estimateof fiber orientation. Accumulated uncertainties in fiber orientation maylead to erroneous reconstruction of pathways. It is difficult to tracereliably in the uncertain areas.

It is desirable to determine and display the tract dispersion, e.g., theeigenvectors and the associated uncertainties. For example, the uniteigenvector may be displayed with a cone of uncertainty around its tip.This conveys the notion that the direction of fiber is not knownprecisely.

However, the methods known in the art are directed to computation andvisualization of a circular cone of uncertainty. These methods are notsuitable for practical computation and visualization of an ellipticalcone of uncertainty.

There is a need for the apparatuses and methods to overcome theabove-referenced problems and others.

BRIEF DESCRIPTION

One embodiment includes a magnetic resonance imaging apparatus,including: a magnetic resonance imaging scanner to acquirediffusion-weighted imaging data of a subject; a reconstruction engine toreconstruct the acquired diffusion-weighted imaging data intodiffusion-weighted image representations of the subject; a diffusiontensor engine to construct a diffusion tensor map of an area of interestof the subject; an eigenvalue/eigenvector ordering engine to obtain andorder pairs of eigenvectors and eigenvalues for at least one voxel inthe diffusion tensor map; a covariance matrix determining engine toconstruct a covariance matrix of a major eigenvector for at least onevoxel for which pairs of eigenvectors and eigenvalues have been obtainedand ordered; a first normalized measure determining engine to compute afirst normalized measure of at least one voxel for which the covariancematrix has been constructed; a second normalized measure determiningengine to compute a second normalized measure of at least one voxel forwhich the covariance matrix has been constructed; and a rendering engineto generate a human-viewable display of an image representation ofuncertainty.

One embodiment includes a magnetic resonance imaging method, including:acquiring a three-dimensional diffusion tensor map of an area ofinterest of a subject; processing the diffusion tensor for at least onevoxel to obtain pairs of eigenvectors and eigenvalues; constructing acovariance matrix of a major eigenvector for at least one voxel;computing first and second normalized uncertainty measures of at leastone voxel; and generating a human-viewable display of an imagerepresentation of uncertainty.

One embodiment includes a magnetic resonance imaging apparatus,including: means for acquiring diffusion-weighted imaging data of asubject; means for reconstructing the acquired diffusion-weightedimaging data into diffusion-weighted image representations of thesubject; means for constructing a diffusion tensor map of an area ofinterest of the subject; means for obtaining and ordering eigenvectorsand eigenvalues for at least one voxel; means for constructing acovariance matrix of a major eigenvector for at least one voxel; meansfor computing a first normalized uncertainty measure; means forcomputing a second normalized uncertainty measure; and means forgenerating a human-viewable display of an image representation ofuncertainty.

One embodiment provides a computer-readable medium comprising software,which software, when executed by a computer system, causes the computerto perform operations to analyze diffusion-weighted imagingrepresentations reconstructed from acquired diffusion-weighted imagingdata of a subject, the computer-readable medium includes: instructionsto construct a diffusion tensor map of an area of interest of thesubject; instructions to obtain and order eigenvectors and eigenvaluesfor at least one voxel; instructions to construct a covariance matrix ofa major eigenvector for at least one voxel; instructions to compute afirst normalized uncertainty measure of at least one voxel; instructionsto compute a second normalized uncertainty measure of at least onevoxel; and instructions to generate a human-viewable display of an imagerepresentation of uncertainty.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other features of various embodiments of the inventionwill be apparent from the following, more particular description of suchembodiments of the invention, as illustrated in the accompanyingdrawings, wherein like reference numbers generally indicate identical,functionally similar, and/or structurally similar elements.

FIG. 1A diagrammatically illustrates an exemplary magnetic resonanceimaging system according to an exemplary embodiment of the invention;

FIG. 1B diagrammatically illustrates a portion of an exemplary magneticresonance imaging system according to an exemplary embodiment of theinvention;

FIG. 1C diagrammatically illustrates the eigenvectors and eigenvalues ofthe diffusion tensor and their relationship with axonal fibers or fiberbundles;

FIG. 2 illustrates a Gnomonic projection according to an exemplaryembodiment of the invention;

FIG. 3 illustrates an inverse Gnomonic projection according to anexemplary embodiment of the invention;

FIG. 4 illustrates an exemplary a cone of uncertainty constructedaccording to an exemplary embodiment of the invention;

FIG. 5A illustrates an exemplary normalized areal measure map;

FIG. 5B illustrates an exemplary normalized circumference map;

FIG. 5C illustrates an exemplary eccentricity map of the ellipse of thecone of uncertainty;

FIG. 5D illustrates an exemplary cones of uncertainty of an axial slicecorresponding to the normalized areal map of FIG. 5A.

DETAILED DESCRIPTION

Exemplary embodiments are discussed in detail below. While specificexemplary embodiments are discussed, it should be understood that thisis done for illustration purposes only. In describing and illustratingthe exemplary embodiments, specific terminology is employed for the sakeof clarity. However, the invention is not intended to be limited to thespecific terminology so selected. A person skilled in the relevant artwill recognize that other components and configurations may be usedwithout parting from the spirit and scope of the invention. It is to beunderstood that each specific element includes all technical equivalentsthat operate in a similar manner to accomplish a similar purpose. Eachreference cited herein is incorporated by reference. The examples andembodiments described herein are non-limiting examples.

With reference to FIG. 1A, a magnetic resonance imaging scanner orsystem 100 may include a housing 102 defining a generally cylindricalscanner bore 104 defining an examination region 105, inside of which anassociated imaging subject 106 is disposed. The subject 106 may be, forexample, a human patient, an animal, a phantom, etc. In FIG. 1A, thehousing 102 is shown in cross-section to illustrate the inside of thehousing 102. Main magnetic field coils 110 may be disposed inside thehousing 102, and may produce a main B₀ magnetic field parallel to acentral axis 112 of the scanner bore 104. In FIG. 1A, the direction ofthe main B₀ magnetic field is parallel to the z-axis of the referencex-y-z Cartesian coordinate system. Main magnetic field coils 110 aretypically superconducting coils disposed inside cryoshrouding 114,although resistive main magnets may also be used. The main magneticfield coils 110 may generate the main B₀ magnetic field, which may besubstantially uniform in an imaging volume of the bore 104.

The housing 102 also houses or supports magnetic field gradient coil(s)116 for selectively producing known magnetic field gradients parallel tothe central axis 112 of the bore 104, along in-plane directionstransverse to the central axis 112, or along other selected directions.In one embodiment, the gradient coil(s) 116 are shielded with shieldingcoil(s) (not shown). The shielding coils are designed to cooperate withthe gradient coil 116 to generate a magnetic field which has asubstantially zero magnetic flux density outside an area defined by theouter radius of the shielding coil(s).

The magnetic resonance imaging scanner 100 may include a radio frequencycoil arrangement or system 120 disposed inside the bore 104 toselectively excite and/or detect magnetic resonances. In one embodiment,the coil system 120 may include a whole body coil 122 and a local coilpositioned in a proximity to the intended volume of examination, such asa head coil 124. The coils 122, 124 may include a TEM coil, a hybridTEM-birdcage coil, a birdcage resonator, an arrangement of loopresonators, or the like. The coils 122, 124 may include a plurality ofresonators positioned around or in the intended volume of examination.The coils 122, 124 may be, for example, circularly cylindrical, but, ofcourse, might have other geometries, such as an elliptic cross-section,semi-circular cross-section, semi-elliptical cross-section, and thelike.

An MRI controller 140 may operate magnetic field gradient controller orcontrollers 142 coupled to the gradients coils 116 to superimposeselected magnetic field gradients on the main magnetic field in theexamination region, and may also operate a radio frequency transmitteror transmitters 144 each coupled to an individual radio frequency coilsegment to inject selected radio frequency excitation pulses at aboutthe magnetic resonance frequency into the examination region forimaging. The radio frequency transmitter or transmitters 144 may beindividually controlled and may have different phases and amplitudes.

Magnetic resonance is generated and spatially encoded in at least aportion of a region of interest of the imaging subject 106. By applyingselected magnetic field gradients via the gradient coils 116, a selectedk-space trajectory is traversed, such as a Cartesian trajectory, aplurality of radial trajectories, or a spiral trajectory. Alternatively,imaging data may be acquired as projections along selected magneticfield gradient directions. During imaging data acquisition, a radiofrequency receiver or receivers 146 coupled to the head coil 124 mayacquire magnetic resonance samples that are stored in a magneticresonance data memory 150. Of course, it is contemplated that themagnetic resonance signals may be excited and received by a single coilarray such as the whole body coil 122 or the local coil 124.

The MRI sequence may include a complex series of magnetic field gradientpulses and/or sweeps generated by the gradient coils 116 which alongwith selected RF pulses generated by the coils 122, 124 result onmagnetic resonance echoes. For diffusion tensor magnetic resonanceimaging (DT-MRI), data is acquired without diffusion weighting, and withdiffusion weighting in N directions. Six diffusion directions mayprovide sufficient information to construct the diffusion tensor at eachvoxel. A reconstruction engine 152 may reconstruct the imaging data intoan image representation, e.g., diffusion-weighted image representations.The reconstructed image generated by the reconstruction engine 152 maybe stored in an image memory 154.

A diffusion tensor engine 160 may calculate a plurality of diffusioncoefficients values on a per voxel basis, as known in the art, toconstruct a diffusion tensor map 162. For each voxel, a tensor, e.g., asymmetric positive definite 3×3 matrix, that describes a threedimensional shape of diffusion may be calculated. The diffusion tensorat each voxel may be processed to obtain eigenvectors and eigenvalues.

As described in a greater detail below, a covariance matrix determiningengine 164 may construct a covariance matrix of, for example, the majoreigenvector at each voxel. It is contemplated that covariance matricesof the medium and minor eigenvectors may be constructed as well. A coneof uncertainty (COU) determining engine 170 may construct or determine acone of uncertainty (COU) based on the covariance matrix. A firstnormalized measure engine 172 may determine a first normalized measureof the elliptical COU. A second normalized measure engine 174 maydetermine a second normalized measure or a normalized circumference ofthe elliptical COU.

A rendering engine 176 may format the constructed cones of uncertainty,normalized areal maps, normalized circumference maps, or any otherappropriate image representations to be displayed on a user interface180, stored in non-volatile memory, transmitted over a local intranet orthe Internet, viewed, stored, manipulated, or so forth.

Major eigenvector maps represent the fiber orientation and may bevisualized as vector-field or color coded maps giving a cartography ofthe tracts' position and direction. The brightness may be weighted by afractional anisotropy which is a scalar measure of a degree ofanisotropy for a given voxel. To visualize the fiber direction map inthe context of a conventional structural image, the fiber map may beregistered to and superimposed on a structural image, for example, on ahigh resolution MR image, as known in the art. The fiber orientationinformation provided by the diffusion imaging maps may be used todifferentiate among nerve fiber pathways to determine abnormalities.Further, the fiber orientation information provided by the diffusionimaging maps may be used to delineate anatomical structures based on apriori distinct fiber architecture of each anatomical structure.

The user interface 180 may also enable a radiologist, technician, orother operator of the magnetic resonance imaging scanner 100 tocommunicate with the magnetic resonance imaging controller 140 toselect, modify, and execute magnetic resonance imaging sequences.

With continuing reference to FIG. 1A and further reference to FIGS. 1Band 1C, an eigenvector/eigenvalue ordering engine 182 may order thediffusion tensor eigenvectors and eigenvalues for each voxel to obtainordered eigenvectors and eigenvalues. For example, theeigenvector/eigenvalue ordering engine 182 may order eigenvalues fromlargest to smallest eigenvalue λ1, λ2, λ3. The eigenvector whichcorresponds to the largest eigenvalue λ1 is called the majoreigenvector, and aligns with the spatial direction having the highestdiffusion coefficient. The medium and smallest eigenvalues λ2, λ3correspond to the medium and minor eigenvectors. The medium and minoreigenvectors are orthogonal to the major eigenvector and align withspatial directions having lower diffusion coefficients. The relativevalues of the eigenvalues λ1, λ2, λ3 are indicative of the spatialorientation and magnitude of the diffusion tensor anisotropy.

As shown in FIG. 1C, the eigenvectors and eigenvalues are geometricallyrepresentable by an ellipsoid 183 whose long axis aligns with the majoreigenvector e1, i.e. with the direction of the highest apparentdiffusion coefficient. The deviation of the ellipsoid 183 from a perfectsphere is representative of the anisotropy of the diffusion tensor. Ananisotropic diffusion coefficient tensor may reflect the influence ofneural fiber bundles 184 which tend to inhibit diffusion in directionspartially or totally orthogonal to the fibers 184, e.g. the directionsof eigenvectors e2, e3. In contrast, diffusion parallel to the fibers184, i.e. channeled along the direction of the major eigenvector e1, isenhanced and larger than along the e2, e3 directions.

Generally, the diffusion tensor, D, is a 3 by 3 symmetric positivedefinite matrix given by:

$\begin{matrix}{D = \begin{pmatrix}D_{xx} & D_{xy} & D_{xz} \\D_{xy} & D_{yy} & D_{yz} \\D_{xz} & D_{yz} & D_{zz}\end{pmatrix}} & \left( {1A} \right)\end{matrix}$

In the discussion below, the measured (noisy) diffusion-weighted signalsare denoted by s_(i). The diffusion-weighted function evaluated at γ,is:

$\begin{matrix}{{{\hat{s}}_{i}(\gamma)} = {\exp\left\lbrack {\sum\limits_{j = 1}^{7}{W_{ij}\gamma_{j}}} \right\rbrack}} & \left( {1B} \right)\end{matrix}$

The diffusion-weighted function evaluated at γ(ξ) is:

$\begin{matrix}{{{\hat{s}}_{i}\left( {\gamma(\xi)} \right)} = {\exp\left\lbrack {\sum\limits_{j = 1}^{7}{W_{ij}{\gamma_{j}(\xi)}}} \right\rbrack}} & \left( {1C} \right)\end{matrix}$where γ is the ordinary representation of the diffusion tensor togetherwith the logarithm of the reference signal:

$\begin{matrix}{{\gamma = {\begin{bmatrix}\gamma_{1} & \gamma_{2} & \gamma_{3} & \gamma_{4} & \gamma_{5} & \gamma_{6} & \gamma_{7}\end{bmatrix}^{T}\mspace{14mu} = \begin{bmatrix}{\ln\left( s_{0} \right)} & D_{xx} & D_{yy} & D_{zz} & D_{xy} & D_{yz} & D_{xz}\end{bmatrix}^{T}}},} & \left( {1D} \right)\end{matrix}$and γ(ξ) is the ordinary representation of the diffusion tensor in termsof the Euler diffusion tensor representation, where the Eulerrepresentation of the diffusion tensor is:

$\begin{matrix}{\xi = {\begin{bmatrix}\xi_{1} & \xi_{2} & \xi_{3} & \xi_{4} & \xi_{5} & \xi_{6} & \xi_{7}\end{bmatrix}^{T}\mspace{14mu} = \begin{bmatrix}{\ln\left( s_{0} \right)} & \lambda_{1} & \lambda_{2} & \lambda_{3} & \theta & \phi & \psi\end{bmatrix}^{T}}} & \left( {1E} \right)\end{matrix}$where s₀ is the parameter for the non-diffusion weighted signal,

-   λ₁, λ₂, and λ₃ denote major, medium, and minor eigenvalues,    respectively, and-   θ, φ, and ψ denote Euler angles.

Suppose that {circumflex over (γ)}_(NLS) is the nonlinear least squareestimate of the parameter vector γ, then the residual vector is a vectorwhose individual component is the difference between the observed andthe expected or estimated signals:r({circumflex over (γ)}_(NLS))=[r ₁({circumflex over (γ)}_(NLS)) . . . r_(n)({circumflex over (γ)}_(NLS))]^(T),  (1F)where r_(i)({circumflex over (γ)}_(NLS))=s_(i)−ŝ_(i)({circumflex over(γ)}_(NLS)) for i=1, . . . ,n.

The design matrix, W, is given by:

$\quad\begin{matrix}\left( \begin{matrix}1 & {{- b_{1}}g_{1x}^{2}} & {{- b_{1}}g_{1y}^{2}} & {{- b_{1}}g_{1z}^{2}} & {{- 2}b_{1}g_{1x}g_{1y}} & {{- 2}b_{1}g_{1y}g_{1z}} & {{- 2}b_{1}g_{1x}g_{1z}} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\1 & {{- b_{n}}g_{nx}^{2}} & {{- b_{n}}g_{ny}^{2}} & {{- b_{n}}g_{nz}^{2}} & {{- 2}b_{n}g_{nx}g_{ny}} & {{- 2}b_{n}g_{ny}g_{nz}} & {{- 2}b_{n}g_{nx}g_{nz}}\end{matrix} \right) & \left( {1G} \right)\end{matrix}$

The objective functions for the nonlinear least squares problem in thediffusion tensor imaging (DTI) in the ordinary diffusion tensorrepresentation and Euler diffusion tensor representation may beexpressed, respectively:

$\begin{matrix}{{{f_{NLS}(\gamma)} = {\frac{1}{2}{\sum\limits_{i = 1}^{n}\underset{\underset{r_{i}{(\gamma)}}{︸}}{\left( {s_{i} - {\exp\left\lbrack {\sum\limits_{j = 1}^{7}{W_{ij}\gamma_{j}}} \right\rbrack}} \right)^{2}}}}},{and}} & \left( {1H} \right) \\{{{f_{ENLS}\left( {\gamma(\xi)} \right)} = {\frac{1}{2}{\sum\limits_{i = 1}^{n}\underset{\underset{r_{i}{({\gamma{(\xi)}})}}{︸}}{\left( {s_{i} - {\exp\left\lbrack {\sum\limits_{j = 1}^{7}{W_{ij}{\gamma_{j}(\xi)}}} \right\rbrack}} \right)^{2}}}}},} & (2)\end{matrix}$

The covariance matrix Σ_(γ) of the ordinary diffusion tensorrepresentation may be defined as:Σ_(γ)=σ_(DW) ² [W ^(T)(Ŝ ² −RŜ)W] ⁻¹,  (3)

The covariance matrix of the Euler diffusion tensor representation maybe derived from the covariance matrix of the ordinary diffusion tensorrepresentation as described below. Generally, the covariance matrixΣ_(ξ) of the Euler diffusion tensor representation may be defined as:Σ_(ξ)=σ_(DW) ² [J _(ξ) ^(T)(γ({circumflex over (ξ)}))W ^(T)(Ŝ ² =RŜ)WJ_(ξ)(γ({circumflex over (ξ)}))]⁻¹  (4A)where the matrix or vector transposition is denoted by superscript T,

-   J_(ξ) denotes a Jacobian matrix,-   S and Ŝ are diagonal matrices including diagonal elements,-   R is a residual matrix, and-   σ_(DW) ² is the estimated variance.

A first diagonal matrix S includes observed diffusion weighted signals:

$\begin{matrix}{S = \begin{pmatrix}s_{1} & \; & \; \\\; & \ddots & \; \\\; & \; & s_{n}\end{pmatrix}} & \left( {4B} \right)\end{matrix}$

A second diagonal matrix Ŝ includes the estimated diffusion weightedsignals, respectively:

$\begin{matrix}{\hat{S} = \begin{pmatrix}{\hat{s}}_{1} & \; & \; \\\; & \ddots & \; \\\; & \; & {\hat{s}}_{n}\end{pmatrix}} & \left( {4C} \right)\end{matrix}$

The residual matrix R is equal to a difference between the first andsecond diagonal matrices:R=S−Ŝ  (4D)

The element of the Jacobian matrix, J_(ξ)(γ) may be defined as:

$\begin{matrix}{{\left\lbrack {J_{\xi}(\gamma)} \right\rbrack_{ij} \equiv \frac{\partial\gamma_{i}}{\partial\xi_{j}}};\mspace{14mu}{and}} & \left( {4E} \right)\end{matrix}$

The estimated variance σ_(DW) ² may be derived from the nonlinear fit(Equation 1H):σ_(DW) ²=2f _(NLS)({circumflex over (γ)}_(NLS))/(n−7).  4F)

The core idea of error propagation is to transform one covariance matrixto another covariance matrix of interest through an appropriate mappingbetween the underlying representations. In DT-MRI, the most fundamentalcovariance matrix is Σ_(γ), from which other covariance matrices may beobtained.

For example, the Euler covariance matrix Σ_(ξ) may be constructed fromthe ordinary covariance matrix Σ_(γ) as follows:Σ_(ξ) =J _(γ)(ξ({circumflex over (γ)}))Σ_(γ) J _(γ) ^(T)(ξ({circumflexover (γ)})),  (5)where J_(γ)(ξ({circumflex over (γ)})),

$\left. {\left\lbrack {J_{\xi}\left( {\xi\left( \hat{\gamma} \right)} \right)} \right\rbrack_{ij} \equiv \frac{\partial{\xi_{i}(\gamma)}}{\partial\gamma_{j}}} \right|_{\gamma = \hat{\gamma}},$is the Jacobian matrix of ξ with respect to γ evaluated at {circumflexover (γ)}.

The Jacobian matrix is a locally linear map from dγ to dξ at {circumflexover (γ)}:dξ=J _(γ)(ξ({circumflex over (γ)}))dγ  (6A)

As shown, the Jacobian matrix between representations plays a criticalrole in transforming one covariance matrix to another and the Jacobianmatrix is, in turn, dependent upon the mapping between representations.

Equations (5) and (4A) are equivalent in principle. For example,equation (4A) may be derived by substituting equation (3) into equation(5):

Σ_(ξ) = J_(γ)(ξ(γ̂))Σ_(γ)J_(γ)^(T)(ξ(γ̂))   = σ_(DW)²J_(γ)(ξ(γ̂))[W^(T)(Ŝ² − RŜ)W]⁻¹J_(γ)^(T)(ξ(γ̂))   = σ_(DW)²(J_(γ)(ξ(γ̂)))⁽⁻¹⁾⁽⁻¹⁾[W^(T)(Ŝ² − RŜ)W]⁻¹(J_(γ)^(T)(ξ(γ̂)))⁽⁻¹⁾⁽⁻¹⁾   = σ_(DW)²[(J_(γ)^(T)(ξ(γ̂)))⁻¹(W^(T)(Ŝ² − RŜ)W)(J_(γ)(ξ(γ̂)))⁻¹]⁻¹   = σ_(DW)²[J_(ξ)^(T)(γ(ξ̂))W^(T)(Ŝ² − RŜ)WJ_(ξ)(γ(ξ̂))]⁻¹

In the above derivation, the following well-known identity is used:J _(ξ)(γ({circumflex over (ξ)}))·J _(γ)(ξ({circumflex over (γ)}))=J_(ξ)(ξ({circumflex over (ξ)}))=I  (6B)

Therefore:[J _(γ)(ξ({circumflex over (γ)}))]⁻¹ =J _(ξ)(γ({circumflex over(ξ)}))  (6C)

In practice, it is easier to compute J_(ξ)(γ({circumflex over (ξ)}))than J_(γ)(ξ({circumflex over (γ)})). An explicit mapping betweenrepresentations is usually needed to construct the Jacobian matrix.Occasionally, such a mapping may not exist or may be too complicated toshed light on the problem at hand. If this is the case then the firstorder matrix perturbation method may be helpful in finding the Jacobianmatrix in the form of dξ=J_(γ)(ξ({circumflex over (γ)}))dγ withoutrequiring an explicit mapping ξ(γ).

The covariance matrix determining engine 164 may construct an ordinarycovariance matrix 185 of the major eigenvector based on the nonlinearleast squares method by using a reformulated first order matrixperturbation method. A notation, a(·,·) is used below to simplify thediscussion.

Let

$\begin{matrix}{{Q = \begin{pmatrix}Q_{11} & Q_{12} & Q_{13} \\Q_{21} & Q_{22} & Q_{23} \\Q_{31} & Q_{32} & Q_{33}\end{pmatrix}},} & \left( {7A} \right) \\{{q_{1} = {\begin{pmatrix}q_{1x} \\q_{1y} \\q_{1z}\end{pmatrix} = \begin{pmatrix}Q_{11} \\Q_{21} \\Q_{31}\end{pmatrix}}},} & \left( {7B} \right) \\{{q_{2} = {\begin{pmatrix}q_{2x} \\q_{2y} \\q_{2z}\end{pmatrix} = \begin{pmatrix}Q_{12} \\Q_{22} \\Q_{32}\end{pmatrix}}},{and}} & \left( {7C} \right) \\{{q_{3} = {\begin{pmatrix}q_{3x} \\q_{3y} \\q_{3z}\end{pmatrix} = \begin{pmatrix}Q_{13} \\Q_{23} \\Q_{33}\end{pmatrix}}},} & \left( {7D} \right)\end{matrix}$where q₁, q₂ and q₃ are major, medium and minor eigenvectors,respectively.

The quadratic form, q_(i) ^(T)·D·q_(j), as a dot product between twovectors, may be presented by:q _(i) ^(T) ·D·q _(j) =a(q _(i) ,q _(j))^(T)·β  (7E)wherea(q _(i) ,q _(j))^(T) ≡[q _(ix) q _(jx) q _(iy) q _(jy) q _(iz) q _(jz)q _(ix) q _(jy) +q _(iy) q _(jx) q _(iy) q _(jz) +q _(iz) q _(jy) q_(ix) q _(jz) +q _(iz) q _(jx)]  (8)andβ≡[D _(xx) D _(yy) D _(zz) D _(xy) D _(yz) D _(xz)].  (9)

The first order matrix perturbation method begins, respectively, withthe eigenvalue equation and the orthonormality condition:D·q _(i)=λ_(i) q _(i),  (10)andq_(i) ^(T)q_(j)=δ_(ij),  (11)where δ_(ij) is the Kronecker delta function, which assumes unity if i=jor is equal to zero otherwise.

Assuming a small variation on both sides of equation (10), e.g.δ(D·q_(i))=δ(λ_(i)q_(i)):δD·q _(i) +D·δq _(i)=δλ_(i) q _(i)+λ_(i) δq _(i) for i, j=1,2,3  (12)

Assuming a small variation on both sides of equation (11), e.g. δ(q_(i)^(T)q_(j))=δ(δ_(ij));δq _(i) ^(T) ·q _(j) +q _(i) ^(T) ·δq _(j)=δ(δ_(ij))=0 for i,j=1,2,3  (13)

Equation (13) implies that:δq _(i) ^(T) ·q _(j) =−q _(i) ^(T) ·δq _(j) for i≠j,  (14)

If q_(i) ^(T)q_(i)=1:q _(i) ^(T) ·δq _(i)=0 or δq _(i) ^(T) ·q _(i)=0 for i=1,2,3  (15)

Taking the dot product between equation (12) with q_(j) ^(T):

$\begin{matrix}{{{{{q_{j}^{T} \cdot \delta}\;{D \cdot q_{i}}} + {{q_{j}^{T} \cdot D \cdot \delta}\; q_{i}}} = {{\delta\;\lambda_{i}\underset{\underset{0}{︸}}{q_{j}^{T} \cdot q_{i}}} + {\lambda_{i}{q_{j}^{T} \cdot \delta}\; q_{i}}}},{{{{q_{j}^{T} \cdot \delta}\;{D \cdot q_{i}}} + {{\left( {D^{T} \cdot q_{j}} \right)^{T} \cdot \delta}\; q_{i}}} = {\lambda_{i}{q_{j}^{T} \cdot \delta}\;{q_{i}.}}}} & (16)\end{matrix}$

Since the diffusion tensor D is symmetric, equation (16) may be furtherreduced by taking into account equation (10):

$\begin{matrix}{{{{{q_{j}^{T} \cdot \delta}\;{D \cdot q_{i}}} + {{\lambda_{j} \cdot q_{j}^{T} \cdot \delta}\; q_{i}}} = {\lambda_{i}{q_{j}^{T} \cdot \delta}\; q_{i}}},} & (17) \\{{{{q_{j}^{T} \cdot \delta}\; q_{i}} = {\frac{1}{\left( {\lambda_{i} - \lambda_{j}} \right)}{q_{j}^{T} \cdot \delta}\;{D \cdot q_{i}}}},{{{q_{j}^{T} \cdot \delta}\; q_{i}} = {\frac{1}{\left( {\lambda_{i} - \lambda_{j}} \right)}{{a\left( {q_{j},q_{i}} \right)}^{T} \cdot \delta}\;\beta}}} & (18)\end{matrix}$

Although equation (18) differs from equation (6A) defining thetransformation relationship between γ and ζ via the Jacobian matrix,equation (18) may be reformulated to achieve the form of equation (6),e.g. δq_(i)=J_(β)(q_(i))δβ.

Assuming that the relationship between the eigenvalues is λ₁>λ₂>λ₃, andsolving for the dispersion of the major eigenvector q₁, e.g., i=1,equation (18) results in equations (19)-(21). Equation (19) may bederived from equation (15):q ₁ ^(T) ·δq ₁=[0 0 0 0 0 0]·δβ  (19)

Equations (20)-(21) may be derived from equations (17)-(18):

$\begin{matrix}{{{{q_{2}^{T} \cdot \delta}\; q_{1}} = {\frac{1}{\left( {\lambda_{1} - \lambda_{2}} \right)}{{a\left( {q_{2},q_{1}} \right)}^{T} \cdot \delta}\;\beta}},} & (20) \\{{{q_{3}^{T} \cdot \delta}\; q_{1}} = {\frac{1}{\left( {\lambda_{1} - \lambda_{3}} \right)}{{a\left( {q_{3},q_{1}} \right)}^{T} \cdot \delta}\;{\beta.}}} & (21)\end{matrix}$

Equations (19-21) may be combined into a single matrix expression:

$\begin{matrix}{{{{{\underset{\underset{Q^{T}}{︸}}{\begin{pmatrix}q_{1}^{T} \\q_{2}^{T} \\q_{3}^{T}\end{pmatrix}} \cdot \delta}\; q_{1}} = {{\underset{\underset{S_{1}}{︸}}{\begin{pmatrix}{0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0} \\{\frac{1}{\left( {\lambda_{1} - \lambda_{2}} \right)}{a\left( {q_{2},q_{1}} \right)}^{T}} \\{\frac{1}{\left( {\lambda_{1} - \lambda_{3}} \right)}{a\left( {q_{3},q_{1}} \right)}^{T}}\end{pmatrix}} \cdot \delta}\;\beta}},{or}}{{{{Q^{T} \cdot \delta}\; q_{1}} = {{S_{1} \cdot \delta}\;\beta}},{or}}{{\delta\; q_{1}} = {{Q \cdot S_{1} \cdot \delta}\;{\beta.}}}} & (22)\end{matrix}$

In the above derivation, the orthonormality condition is used, e.g.QQ^(T)=I. The dispersions of the medium and minor eigenvectors may besimilarly derived.

The above reformulation based on β may be extended to γ resulting in:δq ₁ =Q·T ₁·δγ  (23A)where T₁ is a 3×7 matrix given by

$\begin{matrix}{T_{1} = \begin{pmatrix}{0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0\mspace{31mu} 0} \\{0\mspace{31mu}\frac{1}{\left( {\lambda_{1} - \lambda_{2}} \right)}{a\left( {q_{2},q_{1}} \right)}^{T}} \\{0\mspace{31mu}\frac{1}{\left( {\lambda_{1} - \lambda_{3}} \right)}{a\left( {q_{3},q_{1}} \right)}^{T}}\end{pmatrix}} & \left( {23B} \right)\end{matrix}$

In terms of the Jacobian matrix for β:Q·S ₁ =J _(β)(q ₁)  (24)

In terms of the Jacobian matrix for γ:Q·T ₁ =J _(γ)(q ₁)  (25)

Finally, the covariance matrix Σ_(q) ₁ of the major eigenvector q₁obtained via the reformulated first order matrix perturbation method maybe determined as:Σ_(q) ₁ ({circumflex over (γ)})=J _(γ)(q ₁)Σ_(γ)({circumflex over (γ)})J_(γ) ^(T)(q ₁)  (26)

Optionally or alternatively, the covariance matrix Σ_(q) ₁ derived inequation (26) may be derived based on the error propagation method,known in the art. The covariance matrix of the major eigenvector of thediffusion tensor using the Euler diffusion tensor representation is:Σ_(q) ₁ =J _(ξ)(q ₁({circumflex over (ξ)}))Σ_(ξ) J _(ξ) ^(T)(q₁({circumflex over (ξ)})),  (27)

Substituting equation (4A) into equation (27) results in:

$\begin{matrix}\begin{matrix}{\Sigma_{q_{1}} = {\sigma_{DW}^{2}{{J_{\xi}\left( {q_{1}\left( \hat{\xi} \right)} \right)}\left\lbrack {{J_{\xi}^{T}\left( {\gamma\left( \hat{\xi} \right)} \right)}{W^{T}\left( {{\hat{S}}^{2} - {R\;\hat{S}}} \right)}{{WJ}_{\xi}\left( {\gamma\left( \hat{\xi} \right)} \right)}} \right\rbrack}^{- 1}{J_{\xi}^{T}\left( {q_{1}\left( \hat{\xi} \right)} \right)}}} \\{\Sigma_{q_{1}} = {{J_{\xi}\left( {q_{1}\left( \hat{\xi} \right)} \right)}\left\lbrack {J_{\xi}\left( {\gamma\left( \hat{\xi} \right)} \right)} \right\rbrack}^{- 1}} \\{{\underset{\underset{\Sigma_{\gamma}}{︸}}{{\sigma_{DW}^{2}\left\lbrack {{W^{T}\left( {{\hat{S}}^{2} - {R\;\hat{S}}} \right)}W} \right\rbrack}^{- 1}}\left\lbrack {J_{\xi}^{T}\left( {\gamma\left( \hat{\xi} \right)} \right)} \right\rbrack}^{- 1}{J_{\xi}^{T}\left( {q_{1}\left( \hat{\xi} \right)} \right)}} \\{\Sigma_{q_{1}} = {{J_{\xi}\left( {q_{1}\left( \hat{\xi} \right)} \right)}{J_{\gamma}\left( {\xi\left( \hat{\gamma} \right)} \right)}\Sigma_{\gamma}{J_{\gamma}^{T}\left( {\xi\left( \hat{\gamma} \right)} \right)}{J_{\xi}^{T}\left( {q_{1}\left( \hat{\xi} \right)} \right)}}} \\{\Sigma_{q_{1}} = {{J_{\gamma}\left( {q_{1}\left( \hat{\gamma} \right)} \right)}\Sigma_{\gamma}{J_{\gamma}^{T}\left( {q_{1}\left( \hat{\gamma} \right)} \right)}}}\end{matrix} & \left( {28A} \right)\end{matrix}$where J_(ξ)(q₁({circumflex over (ξ)}))J_(γ)(ξ({circumflex over(γ)}))=J_(γ)(q₁({circumflex over (γ)})), and[J _(ξ)(γ({circumflex over (ξ)}))]⁻¹ =J _(γ)(ξ({circumflex over (γ)})).

Jacobian matrix J_(γ)(q₁({circumflex over (γ)})) may be constructed asdescribed above.

In one embodiment, the cone of uncertainty determining engine 170 mayconstruct the elliptical cone of uncertainty as follows. From equation(28A), the 1−α confidence region for the major eigenvector q₁ is anellipsoid given by:(q ₁ −q ₁({circumflex over (γ)}))^(T)Σ_(q) ₁ ⁺({circumflex over (γ)})(q₁ −q ₁({circumflex over (γ)}))≦2

(2,n−7;α)  (28B)where

(2,n−7;α) is the upper α quantile for an

distribution with 2 and n−7 degrees of freedom.

The plus sign on Σ_(q) ₁ ⁺ denotes the matrix pseudoinverse which isneeded because the matrix rank of Σ_(q) ₁ is two rather than three.

Equivalently:(q ₁ −q ₁({circumflex over (γ)}))^(T)(2

(2,n−7;α)Σ_(q) ₁ ({circumflex over (γ)}))⁺(q ₁ −q ₁({circumflex over(γ)}))≦1  (28C)

Because the eigenvectors of the diffusion tensor are orthogonal and theeigenvectors of Σ_(q) ₁ ({circumflex over (γ)}) are the singularvectors, equation (28C) may be simplified:q ₁ ^(T)(2

(2,n−7;α)Σ_(q) ₁ ({circumflex over (γ)}))⁺ q ₁≦1  (28D)

Let the eigenvalue decomposition of Σ₁ ₁ ({circumflex over (γ)}) be:Σ_(q) ₁ ({circumflex over (γ)})=ω₁ c ₁ c ₁ ^(T)+ω₂ c ₂ c ₂ ^(T)+0q₁({circumflex over (γ)})q ₁ ^(T) ({circumflex over (γ)}),  (28E)where the last term is due to Q·T₁=J_(γ)(q₁) and equation (26).

The principal directions of the ellipse of the COU whose areacorresponds to the 1−α joint confidence region for the major eigenvectorq₁ may be found by selecting q₁ to be equal to one of:√{square root over (2

(2,n−7;α)ω₁)}c ₁ or √{square root over (2

(2,n−7;α)ω₂)}c ₂  (28F)

The term

(2,n−7;α) may be found by solving equation (28G):

$\begin{matrix}{{\alpha - {I_{{({n - 7})}/{({n - 7 + {2}})}}\left( {{\left( {n - 7} \right)/2},1} \right)}} = 0} & \left( {28G} \right)\end{matrix}$where I_(x)(a,b) is the incomplete Beta function.

With continuing reference to FIGS. 1A and 1B and further reference toFIG. 2, the first normalized measure determining engine 170 maydetermine the normalized first or areal measure Γ. More specifically, asimple closed or Jordan curve 200 on the unit sphere 210 divides theunit sphere into first and second regions 220, 230. For example, thefirst region 220 is greater than the second region 230. The firstnormalized measure determining engine 172 may determine the normalizedareal measure Γ as a ratio of an area of the second, smaller region 230on the unit sphere, that is enclosed by the simple closed curve 200whose Gnomonic or central projection is on the u, v-plane 240, to anarea of a hemisphere.

Let s be a point on an upper hemisphere whose Gnomonic projection on theu, v-plane 240 is the point p. Let p be denoted in a vector form by:

$\begin{matrix}{p = \begin{pmatrix}u \\v \\1\end{pmatrix}} & (29)\end{matrix}$

The point s in a vector form may be obtained by normalizing p to a unitlength:

$\begin{matrix}{{s \equiv \begin{pmatrix}X_{s} \\Y_{s} \\Z_{s}\end{pmatrix}} = {\frac{p}{\sqrt{p^{T} \cdot p}} = \begin{pmatrix}\frac{u}{\sqrt{u^{2} + v^{2} + 1}} \\\frac{v}{\sqrt{u^{2} + v^{2} + 1}} \\\frac{1}{\sqrt{u^{2} + v^{2} + 1}}\end{pmatrix}}} & \left( {30A} \right)\end{matrix}$

Therefore, the components of s are functions of u and v.

Suppose that major and minor axes 242, 244 of the ellipse 246 coincidewith the u and the v axes of the u, v-plane. Analogously to equation(28F), let a length of the major and the minor axes 242, 244 be a and b,respectively, where:a=√{square root over (2

(2,n−7;α)ω₁)} and  (30B)b=√{square root over (2

(2,n−7;α)ω₂)}  (30C)

The normalized areal measure, Γ, which is a dimensionless quantity, hasa range of values from zero to unity:

$\begin{matrix}{\Gamma = {\frac{1}{2\pi}{\int{\int_{R}\ {\sqrt{{EG} - F^{2}}{\mathbb{d}u}{\mathbb{d}v}}}}}} & \left( {31A} \right)\end{matrix}$where R is the smaller region 230.

A region R may be is defined by:(u/a)²+(v/b)²≦1  (31B)

The coefficients of the first fundamental form are:

$\begin{matrix}{{E = {\left( \frac{\partial X_{s}}{\partial u} \right)^{2} + \left( \frac{\partial Y_{s}}{\partial u} \right)^{2} + \left( \frac{\partial Z_{s}}{\partial u} \right)^{2}}},} & \left( {31C} \right) \\{{G = {\left( \frac{\partial X_{s}}{\partial v} \right)^{2} + \left( \frac{\partial Y_{s}}{\partial v} \right)^{2} + \left( \frac{\partial Z_{s}}{\partial v} \right)^{2}}},{and}} & \left( {31D} \right) \\{F = {{\frac{\partial X_{s}}{\partial u}\frac{\partial X_{s}}{\partial v}} + {\frac{\partial Y_{s}}{\partial u}\frac{\partial Y_{s}}{\partial v}} + {\frac{\partial Z_{s}}{\partial u}{\frac{\partial Z_{s}}{\partial v}.}}}} & \left( {31E} \right)\end{matrix}$

The surface integral of equation (31A) may be simplified to showexplicit dependence on a and b. First, several preliminary expressionsare evaluated:

$\begin{matrix}{{\frac{\partial X_{s}}{\partial u} = \frac{1 + v^{2}}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}},} & \left( {31F} \right) \\{{\frac{\partial X_{s}}{\partial v} = {\frac{uv}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}} = \frac{\partial Y_{s}}{\partial u}}},} & \left( {31G} \right) \\{{\frac{\partial Y_{s}}{\partial v} = \frac{1 + u^{2}}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}},} & \left( {31H} \right) \\{{\frac{\partial Z_{s}}{\partial u} = {- \frac{u}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}}},} & \left( {31I} \right) \\{{\frac{\partial Z_{s}}{\partial v} = {- \frac{v}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}}};} & \left( {31\; J} \right) \\{{E = {{\left( \frac{\partial X_{s}}{\partial u} \right)^{2} + \left( \frac{\partial Y_{s}}{\partial u} \right)^{2} + \left( \frac{\partial Z_{s}}{\partial u} \right)^{2}} = \frac{1 + v^{2}}{\left( {u^{2} + v^{2} + 1} \right)^{2}}}},} & \left( {31K} \right) \\{{G = {{\left( \frac{\partial X_{s}}{\partial v} \right)^{2} + \left( \frac{\partial Y_{s}}{\partial v} \right)^{2} + \left( \frac{\partial Z_{s}}{\partial v} \right)^{2}} = \frac{1 + u^{2}}{\left( {u^{2} + v^{2} + 1} \right)^{2}}}},{and}} & \left( {31L} \right) \\{F = {{{\frac{\partial X_{s}}{\partial u}\frac{\partial X_{s}}{\partial v}} + \frac{\partial Y_{s}}{\partial u} + {\frac{\partial Z_{s}}{\partial u}\frac{\partial Z_{s}}{\partial v}}} = {- {\frac{uv}{\left( {u^{2} + v^{2} + 1} \right)^{2}}.}}}} & \left( {31M} \right)\end{matrix}$

The integrand of equation (31A) is then given by:

$\begin{matrix}{\sqrt{{EG} - F^{2}} = {\sqrt{\frac{\left( {1 + v^{2}} \right)\left( {1 + u^{2}} \right)}{\left( {u^{2} + v^{2} + 1} \right)^{4}} - \frac{u^{2}v^{2}}{\left( {u^{2} + v^{2} + 1} \right)^{4}}}\mspace{121mu} = \frac{1}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}}} & \left( {31N} \right)\end{matrix}$

Therefore,

$\begin{matrix}{\Gamma = {\frac{1}{2\pi}{\int{\int_{R}\ {\frac{1}{\left( {u^{2} + v^{2} + 1} \right)^{3/2}}\ {\mathbb{d}u}{\mathbb{d}v}}}}}} & \left( {31O} \right)\end{matrix}$where R is the region defined by (u/a)²+(v/b)²≦1.

By a change of variables, u=aξ_(x) and v=bξ_(x):

$\begin{matrix}{\Gamma = {\frac{ab}{2\;\pi}{\int{\int_{R_{1}}{\frac{1}{\left( {1 + {a^{2}\xi_{x}^{2}} + {b^{2}\xi_{y}^{2}}} \right)^{3/2}}\ {\mathbb{d}\xi_{x}}{\mathbb{d}\xi_{y}}}}}}} & \left( {31P} \right)\end{matrix}$where R₁ is a circular region (or a disk) defined by ξ_(x) ²+ξ_(y) ²≦1

The limits of integration due to the circular region R₁ may beintroduced into Γ and the integral is now given by:

$\begin{matrix}{\Gamma = {\frac{ab}{2\;\pi}{\int_{- 1}^{1}{\int_{- \sqrt{1 - \xi_{y}^{2}}}^{\sqrt{1 - \xi_{y}^{2}}}{\frac{1}{\left( {1 + {a^{2}\xi_{x}^{2}} + {b^{2}\xi_{y}^{2}}} \right)^{3/2}}\ {\mathbb{d}\xi_{x}}\ {\mathbb{d}\xi_{y}}}}}}} & \left( {31Q} \right)\end{matrix}$

In equation (31Q), the term a² may be factored and β² may be defined as:

$\begin{matrix}{\beta^{2} \equiv {\frac{1}{a^{2}} + {\frac{b^{2}}{a^{2}}\xi_{y}^{2}}}} & \left( {31R} \right)\end{matrix}$

The inner integral of equation (31Q) becomes:

$\begin{matrix}{\int{\frac{1}{\left( {\beta^{2} + \xi^{2}} \right)^{3/2}}{\mathbb{d}\xi}}} & \left( {31S} \right)\end{matrix}$

The integral of equation (31S) may be evaluated as described below.Define

$\begin{matrix}{T = {\xi\left( {\beta^{2} + \xi^{2}} \right)}^{{- \frac{m}{2}} + 1}} & \left( {31T} \right)\end{matrix}$

The derivative of T with respect to ξ is given by:

$\begin{matrix}\begin{matrix}{\frac{\mathbb{d}T}{\mathbb{d}\xi} = {\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1} - {{\xi^{2}\left( {m - 2} \right)}\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}}}}} \\{= {\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1} - {\left( {\left( {\beta^{2} + \xi^{2}} \right) - \beta^{2}} \right)\left( {m - 2} \right)\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}}}}} \\{= {\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1} - {\left( {m - 2} \right)\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1}} +}} \\{{\beta^{2}\left( {m - 2} \right)}\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}}} \\{= {{\left( {3 - m} \right)\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1}} + {{\beta^{2}\left( {m - 2} \right)}{\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}}.}}}}\end{matrix} & \left( {31U} \right)\end{matrix}$

Rearranging β² (m−2) in equation (31U):

$\begin{matrix}{\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}} = {{\frac{\left( {3 - m} \right)}{\beta^{2}\left( {m - 2} \right)}{\int\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1}}} + {\frac{1}{\beta^{2}\left( {m - 2} \right)}\frac{\mathbb{d}T}{\mathbb{d}\xi}}}} & \left( {31V} \right)\end{matrix}$

Therefore,

$\begin{matrix}{{{\int{\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{m}{2}}{\mathbb{d}\xi}}} = {{\frac{\left( {3 - m} \right)}{\beta^{2}\left( {m - 2} \right)}{\int{\left( {\beta^{2} + \xi^{2}} \right)^{{- \frac{m}{2}} + 1}{\mathbb{d}\xi}}}} + {\frac{1}{\beta^{2}\left( {m - 2} \right)}{\xi\left( {\beta^{2} + \xi^{2}} \right)}^{{- \frac{m}{2}} + 1}}}}\mspace{85mu}{or}\mspace{79mu}{{\int{\left( {\beta^{2} + \xi^{2}} \right)^{- \frac{3}{2}}{\mathbb{d}\xi}}} = {\frac{1}{\beta^{2}}{{\xi\left( {\beta^{2} + \xi^{2}} \right)}^{\frac{1}{2}}.}}}} & \left( {31\; W} \right)\end{matrix}$

In one embodiment, for a faster approach, a trigonometric substitutionof tan θ=ξ/β may be used.

The areal measure Γ:

$\begin{matrix}{{\Gamma = {\frac{b}{2\;\pi\; a^{2}}{\int_{- 1}^{1}{\left\lbrack {\frac{1}{\beta^{2}}{\xi_{x}\left( {\beta^{2} + \xi_{x}^{2}} \right)}^{- \frac{1}{2}}} \right\rbrack_{\xi_{x} = {- \sqrt{1 - \xi_{y}^{2}}}}^{\xi_{x} = {+ \sqrt{1 - \xi_{y}^{2}}}}\ {\mathbb{d}\xi_{y}}}}}},} \\{{= {\frac{b}{\pi\; a^{2}}{\int_{- 1}^{1}{\left\lbrack {\frac{1}{\beta^{2}}\sqrt{1 - \xi_{y}^{2}}\left( {\beta^{2} + 1 - \xi_{y}^{2}} \right)^{- \frac{1}{2}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}}},} \\{{= {\frac{ab}{\pi}{\int_{- 1}^{1}{\left\lbrack {\frac{\sqrt{1 - \xi_{y}^{2}}}{1 + {b^{2}\xi_{y}^{2}}}\ \frac{1}{\left( {1 + {b^{2}\xi_{y}^{2}} + a^{2} - {a^{2}\xi_{y}^{2}}} \right)^{1/2}}} \right\rbrack{\mathbb{d}\xi_{y}}}}}},}\end{matrix}$

${= {\frac{2\;{ab}}{\pi}{\int_{0}^{1}{\left\lbrack {\frac{\sqrt{1 - \xi_{y}^{2}}}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{\left( {1 + {b^{2}\xi_{y}^{2}} + a^{2} - {a^{2}\xi_{y}^{2}}} \right)^{1/2}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}}},$since the integrand is an even function,

$\begin{matrix}{= {\frac{2\;{ab}}{\pi\sqrt{1 + a^{2}}}{\int_{0}^{1}{\left\lbrack {\frac{1 - \xi_{y}^{2}}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{\begin{matrix}\left( {1 - \xi_{y}^{2}} \right)^{1/2} \\\left( {1 - {\frac{a^{2} - b^{2}}{1 + a^{2}}\xi_{y}^{2}}} \right)^{1/2}\end{matrix}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}}} & \left( {31X} \right)\end{matrix}$

The integral I of equation (31X) may be expressed in terms of thecomplete elliptic integrals of the first and third kinds.

Define:

$\begin{matrix}{{d_{1} = {\frac{1}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{\left( {1 - \xi_{y}^{2}} \right)^{1/2}\left( {1 - {\frac{a^{2} - b^{2}}{1 + a^{2}}\xi_{y}^{2}}} \right)^{1/2}}}}{and}} & \left( {31Y} \right) \\{{{d_{2} = \frac{1}{\left( {1 - \xi_{y}^{2}} \right)^{1/2}\left( {1 - {\frac{a^{2} - b^{2}}{1 + a^{2}}\xi_{y}^{2}}} \right)^{1/2}}};}{{then}\text{:}}} & \left( {31Z} \right) \\{{{K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)} = {\int_{0}^{1}{\frac{1}{d_{2}}\ {\mathbb{d}\xi_{y}}}}}{and}} & \left( {31{AA}} \right) \\{{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} = {\int_{0}^{1}\ {\frac{1}{d_{1}}{\mathbb{d}\xi_{y}}}}} & \left( {31{BB}} \right)\end{matrix}$

Therefore, the integral I may be expressed in two formats as shown inequations (31CC) and (31DD):

$\begin{matrix}{I = {\int_{0}^{1}{\left\lbrack {\frac{1 - \xi_{y}^{2}}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{d_{2}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}} & \left( {31{CC}} \right) \\{I = {\int_{0}^{1}{\left\lbrack \frac{1 - \xi_{y}^{2}}{d_{1}} \right\rbrack{\mathbb{d}\xi_{y}}}}} & \left( {31{DD}} \right)\end{matrix}$

For the integral I of equation (31CC):

$\begin{matrix}\begin{matrix}{{I = {\int_{0}^{1}{\left\lbrack {\frac{1 - \xi_{y}^{2} + {b^{2}\xi_{y}^{2}} - {b^{2}\xi_{y}^{2}}}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{d_{2}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}},} \\{{= {\int_{0}^{1}{\left\lbrack {\frac{1 + {b^{2}\xi_{y}^{2}} - {\left( {1 + b^{2}} \right)\xi_{y}^{2}}}{1 + {b^{2}\xi_{y}^{2}}}\frac{1}{d_{2}}} \right\rbrack\ {\mathbb{d}\xi_{y}}}}},} \\{{= {{\int_{0}^{1}{\frac{1}{d_{2}}\ {\mathbb{d}\xi_{y}}}} - {\left( {1 + b^{2}} \right){\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}}},} \\{= {{K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)} - {\left( {1 + b^{2}} \right){\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}}}\end{matrix} & \left( {31{EE}} \right)\end{matrix}$

For the integral I of equation (31DD):

$\begin{matrix}{I = {{\int_{0}^{1}{\left\lbrack \frac{1 - \xi_{y}^{2}}{d_{1}} \right\rbrack\ {\mathbb{d}\xi_{y}}}} = {{{\int_{0}^{1}{\frac{1}{d_{1}}\ {\mathbb{d}\xi_{y}}}} - {\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}\mspace{14mu} = {{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} - {\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {{\mathbb{d}\xi_{y}}.}}}}}}} & \left( {31{FF}} \right)\end{matrix}$

Since the expressions (31EE) and (31FF) are equivalent, the expressions(31EE) and (31FF) may be equated:

$\begin{matrix}{{{{K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)} - {\left( {1 + b^{2}} \right){\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}} = {{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} - {\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}}\mspace{79mu}{or}} & \left( {31\;{GG}} \right) \\{\mspace{79mu}{{\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}} = {{- \frac{1}{b^{2}}}\left( {{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} - {K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}} \right)}}} & \left( {31{HH}} \right)\end{matrix}$

The integral, I, may be expressed in terms of K and Π:

$\begin{matrix}{I = {{{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} - {\int_{0}^{1}{\frac{\xi_{y}^{2}}{d_{1}}\ {\mathbb{d}\xi_{y}}}}}\mspace{14mu} = {{{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} + \mspace{40mu}{\frac{1}{b^{2}}\left( {{\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)} - {K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}} \right)}}\mspace{14mu} = {{\left( {1 + \frac{1}{b^{2}}} \right){\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)}} - {\frac{1}{b^{2}}{K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}}}}}} & \left( {31\;{II}} \right)\end{matrix}$

The areal measure Γ may be found from equations (31X) and (31II) as:

$\begin{matrix}{{\Gamma = {\frac{2\;{ab}}{\pi\sqrt{1 + a^{2}}}\begin{pmatrix}{{\left( {1 + \frac{1}{b^{2}}} \right){\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)}} -} \\{\frac{1}{b^{2}}{K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}}\end{pmatrix}}}{or}{{\Gamma = {\frac{2a}{\pi\; b\sqrt{1 + a^{2}}}\left\lbrack {{\left( {1 + b^{2}} \right){\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)}} - {K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}} \right\rbrack}},}} & \left( {32\; A} \right)\end{matrix}$where K and Π are, respectively, the complete elliptic integrals of thefirst kind and the third kind.

The elliptic integral K of the first kind may be defined as:

$\begin{matrix}{{K(m)} = {{\int_{0}^{\pi/2}{\frac{1}{\sqrt{1 - {m\;\sin^{2}\theta}}}\ {\mathbb{d}\theta}}}\mspace{56mu} = {\int_{0}^{1}{\frac{1}{\sqrt{\left( {1 - v^{2}} \right)\left( {1 - {mv}^{2}} \right)}}\ {\mathbb{d}v}}}}} & \left( {32B} \right)\end{matrix}$

The elliptic integral Π of the third kind may be defined as:

$\begin{matrix}{{\Pi\left( {n,m} \right)} = {{\int_{0}^{\pi/2}{\frac{1}{\left( {1 - {n\;\sin^{2}\theta}} \right)\sqrt{1 - {m\;\sin^{2}\theta}}}\ {\mathbb{d}\theta}}}\mspace{85mu} = {\int_{0}^{1}{\frac{1}{\left( {1 - {nv}^{2}} \right)\left( {1 - {mv}^{2}} \right)}\ {\mathbb{d}v}}}}} & \left( {32C} \right)\end{matrix}$

In one embodiment, where the ellipse on the u, v-plane is becomes acircle, e.g., r≡a=b, the normalized areal measure is:

$\begin{matrix}{\Gamma = {1 - {\frac{1}{\sqrt{1 + r^{2}}}.}}} & (33)\end{matrix}$

A normalized areal measure map engine 177 may construct a normalizedareal measure map based on the determined normalized areal measure of anarea of interest. The constructed normalized areal measurement map maybe used by the rendering engine 176 to generate a human-viewable displayof an image of uncertainty which may be displayed on the display 180.

With continuing reference to FIGS. 1B and 2 and further reference toFIG. 3, the second normalized measure determining engine 174 maydetermine the normalized second or circumferential measure as a ratio ofthe circumference of the simple closed curve 200 on the unit sphere 210to the circumference of a great circle 260 of the unit sphere 210. Asexplained in detail below, the normalized circumference of a simpleclosed curve on the sphere may be expressed in terms of the completeelliptical integrals of the first and third kinds. More specifically,the normalized circumference of an inverse projected curve 330 on theunit sphere 210 may be determined as:

$\begin{matrix}{\Lambda = {\frac{1}{2\;\pi}{\int_{0}^{2\;\pi}{\sqrt{\frac{\mathbb{d}{X_{s}(\theta)}}{\mathbb{d}(\theta)^{2}} + \left( \frac{\mathbb{d}\;{Y_{s}(\theta)}}{\mathbb{d}\;\theta} \right)^{2} + \left( \frac{\mathbb{d}{Z_{s}(\theta)}}{\mathbb{d}\;\theta} \right)^{2}}{\mathbb{d}\theta}}}}} & \left( {34\; A} \right)\end{matrix}$

The 2π term in the denominator is the normalization factor of thecircumference of the great circle of the unit sphere. The components ofthe vector, s, may be parameterized in terms of a, b, and θ:

$\begin{matrix}{{{s(\theta)} \equiv \begin{pmatrix}{X_{s}(\theta)} \\{Y_{s}(\theta)} \\{Z_{s}(\theta)}\end{pmatrix}} = \begin{pmatrix}\frac{a\;{\cos(\theta)}}{\sqrt{\left( {a\;{\cos(\theta)}} \right)^{2} + \left( {b\;{\sin(\theta)}} \right)^{2} + 1}} \\\frac{b\;{\sin(\theta)}}{\sqrt{\left( {a\;{\cos(\theta)}} \right)^{2} + \left( {b\;{\sin(\theta)}} \right)^{2} + 1}} \\\frac{1}{\sqrt{\left( {a\;{\cos(\theta)}} \right)^{2} + \left( {b\;{\sin(\theta)}} \right)^{2} + 1}}\end{pmatrix}} & \left( {34B} \right)\end{matrix}$

Equation (34A) may be further simplified and expressed in terms of thecomplete elliptic integrals of the first and third kinds. Provided hereare several preliminary expressions:

$\begin{matrix}{\left( \frac{\mathbb{d}{X_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2} = \frac{{a^{2}\left( {1 + b^{2}} \right)}^{2}{\sin^{2}(\theta)}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack^{3}}} & \left( {34C} \right) \\{\left( \frac{\mathbb{d}{Y_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2} = \frac{{b^{2}\left( {1 + a^{2}} \right)}^{2}{\cos^{2}(\theta)}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack^{3}}} & \left( {34D} \right) \\{\left( \frac{\mathbb{d}{Z_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2} = \frac{\left( {b^{2} - a^{2}} \right)^{2}{\sin^{2}(\theta)}{\cos^{2}(\theta)}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack^{3}}} & \left( {34E} \right) \\{\Phi \equiv {\left( \frac{\mathbb{d}{X_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2} + \left( \frac{\mathbb{d}{Y_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2} + \left( \frac{\mathbb{d}{Z_{s}(\theta)}}{\mathbb{d}\theta} \right)^{2}}} & \left( {34F} \right) \\{\Phi = \frac{\begin{matrix}{{{a^{2}\left( {1 + b^{2}} \right)}^{2}{\sin^{2}(\theta)}} + {{b^{2}\left( {1 + a^{2}} \right)}^{2}\cos^{2}(\theta)} +} \\{\left( {b^{2} - a^{2}} \right)^{2}{\sin^{2}(\theta)}{\cos^{2}(\theta)}}\end{matrix}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack^{3}}} & \left( {34E} \right)\end{matrix}$

Equation (34E) may be expressed as:

$\begin{matrix}{\Phi = \frac{\begin{matrix}\left( {{b^{2}\left( {1 + a^{2}} \right)} + {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right) \\\left( {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right)\end{matrix}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack^{3}}} & \left( {34F} \right)\end{matrix}$

The square root of Φ may be written as:

$\begin{matrix}{\mspace{79mu}{{\sqrt{\Phi} = \frac{\sqrt{{b^{2}\left( {1 + a^{2}} \right)} + {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack}}\mspace{79mu}{or}}} & \left( {34F} \right) \\{\sqrt{\Phi} = {\frac{{b^{2}\left( {1 + a^{2}} \right)} + {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}}{\left\lbrack {1 + a^{2} - {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}} \right\rbrack\sqrt{{b^{2}\left( {1 + a^{2}} \right)} + {\left( {a^{2} - b^{2}} \right){\sin^{2}(\theta)}}}}\mspace{45mu} = {\frac{b}{\left( {1 + a^{2}} \right)^{1/2}}{\frac{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}.}}}} & \left( {34G} \right)\end{matrix}$

The normalized circumference is then given by:

$\begin{matrix}{\Lambda = {\frac{1}{2\;\pi}{\int_{0}^{2\pi}{\sqrt{\Phi}\ {\mathbb{d}\theta}}}}} & \left( {34H} \right)\end{matrix}$

Since an arc length for each quadrant is the same, the normalizedcircumference, Λ, may be reduced to:

$\begin{matrix}{\mspace{79mu}{{\Lambda = {\frac{2}{\pi}{\int_{0}^{\pi/2}{\sqrt{\Phi}\ {\mathbb{d}\theta}}}}}\mspace{14mu}\mspace{79mu}{or}}} & \left( {34I} \right) \\{\Lambda = {\frac{2b}{{\pi\left( {1 + a^{2}} \right)}^{1/2}}{\int_{0}^{\pi/2}{\frac{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}}}} & \left( {34J} \right)\end{matrix}$

The integral I₁ of equation (34J) may be expressed in two equivalentequations (34K) and (34N):

$\begin{matrix}{I_{1} = {\int_{0}^{\pi/2}{\frac{\begin{matrix}{1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}} + {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}\sin^{2}(\theta)} -} \\{\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}\end{matrix}}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}}} & \left( {34\; K} \right) \\{\mspace{20mu}{= {{\int_{0}^{\pi/2}{\frac{1}{\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}} + {\left( {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)} - \frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}} \right)\underset{\underset{II}{︸}}{\int_{0}^{\pi/2}{\frac{\sin^{2}(\theta)}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}}}}}} & \left( {34\; L} \right) \\{\mspace{31mu}{= {{K\left( \frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)} \right)} + {\left( {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)} - \frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}} \right){II}}}}} & \left( {34M} \right) \\{{I_{2} = {\int_{0}^{\pi/2}{\frac{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}}},} & \left( {34N} \right) \\{\mspace{31mu}{= {{\int_{0}^{\pi/2}{\frac{1}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}} - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}\underset{\underset{II}{︸}}{\int_{0}^{\pi/2}{\frac{\sin^{2}(\theta)}{\left\lbrack {1 - {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}} \right\rbrack\sqrt{1 - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{\sin^{2}(\theta)}}}}\ {\mathbb{d}\theta}}}}}}} & \left( {34O} \right) \\{\mspace{31mu}{= {{\Pi\left( {\frac{\left( {a^{2} - b^{2}} \right)}{\left( {1 + a^{2}} \right)},\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}} \right)} - {\frac{\left( {b^{2} - a^{2}} \right)}{b^{2}\left( {1 + a^{2}} \right)}{II}}}}} & \left( {34P} \right)\end{matrix}$

Since the two integral expressions, I₁ and I₂, are equivalent, II may besolved by equating the two expressions so that II is given in terms ofa, b, K and Π. Substituting the expression of II from equation (34P)into equation (34M), the normalized circumference is:

$\begin{matrix}{{\Lambda\left( {a,b} \right)} = {\frac{2}{\pi\; b\sqrt{1 + a^{2}}}\begin{pmatrix}{{\left( {1 + b^{2}} \right){\Pi\left( {\frac{a^{2} - b^{2}}{1 + a^{2}},\frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}}} \right)}} -} \\{K\left( \frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}} \right)}\end{pmatrix}}} & \left( {34Q} \right)\end{matrix}$

The normalized circumferential measure is a dimensionless quantity. Inone embodiment, where the ellipse is a circle, e.g., r≡a=b, thenormalized circumferential measure may be reduced to:

$\begin{matrix}{\Lambda = \frac{r}{\sqrt{1 + r^{2}}}} & \left( {34R} \right)\end{matrix}$

A normalized circumference measure map engine 178 may construct anormalized circumference map based on the determined normalizedcircumferential measure of an area of interest. The constructednormalized circumference map may be used by the rendering engine 176 togenerate a human-viewable display of an image of uncertainty which maybe displayed on the display 180.

With continuing reference to FIGS. 1B, and 2-3 and further reference toFIG. 4, the COU determining engine 170 may construct a cone ofuncertainty 410 based on the inverse Gnomonic projection of an ellipse340 of the u, v-plane 240 into the unit sphere 210, illustrated in FIG.3.

In an example provided below, the technique of error propagation viadiffusion tensor representations is denoted by “EP” while thereformulated perturbation technique described above is denoted by “RP”.

A synthetic diffusion tensor, as known in the art, may be given by:γ=[ln(1000)×10⁺⁴ s/mm² 9.475 6.694 4.829 1.123 −0.507 −1.63]^(T)×10⁻⁴mm²/s

The eigenvalue-eigenvector pairs of the example tensor are:{10.4×10⁻⁴ mm²/s,q ₁=[0.9027 0.3139 −0.2940]^(T)},{6.30×10⁻⁴ ,q ₂=[0.3197 −0.9470 −0.0295]^(T)}, and{4.30×10⁻⁴ ,q ₃=[0.2877 0.0673 0.9553]^(T)}.

Further, the trace and fractional anisotropy (FA) are 0.0021 mm²/s and0.4176, respectively.

The mean covariance matrices of the major eigenvector of the exampletensor, based on EP and RP, are substantially equal to one another:

${{{}_{}^{}{}_{q1}^{}} \approx {{}_{}^{}{}_{q1}^{}}} = {\begin{pmatrix}6.0911 & {- 13.269} & 4.5350 \\{- 13.269} & 40.379 & 2.3675 \\4.5350 & 2.3675 & 16.450\end{pmatrix} \times {10^{- 5}.}}$

The computation above is carried out based on a signal-to-noise ratio(SNR) of 50 and on a design matrix, W, that was constructed from a 35gradient direction set with four spherical shells having b values of 0,500, 1000, and 1500 s/mm².

The eigenvalue-eigenvector pairs of _(EP)Σ₁ ₁ are:{ω₁=4.4935×10⁻⁴,[0.3202, −0.9469, −0.0277]^(T)},{ω₂=1.7985×10⁻⁴,[0.2871 0.0691 0.9553]^(T)},and{ω₃=4.387×10⁻²⁰≅0,[0.9027 0.3139 −0.2940]^(T)}.

The only difference between the eigenvalue-eigenvector pairs of_(EP)Σ_(q) ₁ and of _(RP)Σ_(q) ₁ is in the minor eigenvalue,2.5243×10⁻²⁰ for _(RP)Σ_(q) ₁ , which is of no consequence in practice.

The angle of deviation of the two eigenvectors of _(RP)Σ_(q) ₁ may becompared to that of the circular cone of uncertainty. Since the elementsin δD of equation (26) are taken to be the standard deviation of thediffusion tensor elements, let α be equal to approximately 0.3173, whichis the area covering the two tails of the normal distribution at onestandard deviation apart from the center, i.e., at −1 and +1. Therefore:

(2,140−7;0.3173)=1.1578

The angles of deviation for the major and the minor eigenvectors of_(RP)Σ_(q) ₁ are given by θ₁=tan⁻¹(√{square root over (2

(2,n−7;α)ω₁)}), and θ₂=tan⁻¹(√{square root over (2

(2,n−7;α)ω₂)}), respectively, where ω₁ and ω₂ are the first two largesteigenvalues of _(RP)Σ_(q) ₁ . In the context of the example discussedhere, θ₁=1.847° and θ₂=1.169°.

With reference to FIGS. 5A, 5B, 5C and 5D, using experimental simulatedhuman head tensor data together with the experimental design defined inthe above example, the normalized areal measure map 510, whichcorresponds to the 0.95 joint confidence region (or 95% confidenceregion) at an SNR level of 15 may be obtained. Likewise, the normalizedcircumference map 520 may be obtained. The eccentricity map 530 of theellipse of the cone of uncertainty is given by √{square root over(1−b²/a²)}. It is assumed that b≦a. The cones of uncertainty 540 of anaxial slice corresponding to the normalized areal map 510 areillustrated in FIG. 5D.

As described above, the perturbation method is reformulated to obtainthe covariance of the major eigenvector of the diffusion tensor. Thecovariance matrix is shown to be equivalent to that derived from theerror propagation method based on the Euler representation. When amapping between representations of interest is not available then it islikely that the first order matrix perturbation method may be helpful infinding the Jacobian matrix so that the transformation from onecovariance structure to another may be realized. Finally, two normalizedmeasures of the cone of uncertainty and a technique of visualizing coneof uncertainty are described.

The proposed measures, which are directly linked to the uncertainty inthe major eigenvector of the diffusion tensor, may be important forprobing the integrity of the white matter tracts in the brain. Inaddition to that, the measures have a dependence on signal-to-noiseratio and may be used as a calibration gauge for assessing MRI system orDT-MRI acquisition performance. The described measures for quantifyinguncertainty of the major eigenvector of the diffusion tensor aredimensionless and normalized to unity. The measures have directgeometric interpretations. In particular, the areal measure correspondsdirectly to the projection of the 1+α joint confidence region for q₁ inthe u, v-plane onto the unit sphere. Besides the areal measure, thecircumferential measure may be important in gaining insights into theasymmetric nature of tract dispersion. The dispersion informationcontained in Σ_(q) ₁ may be incorporated into the modeling the priordistribution of fiber coherence in probabilistic (or Bayesian)tractography.

By using the techniques described above, the COU may be constructedwithout overlapping cones in neighboring regions.

In an exemplary embodiment, the engines 152, 160, 164, 170, 172, 174,176, 177, 178, and 182 may be implemented with, for example, a computerhaving a processing unit, persistent memory, non-persistent memory,input device(s), output device(s), etc. The computer may includesoftware stored in memory, which when executed by the computer, causesthe computer to perform the functions of one or more of the engines. Thecomputer may include one or more dedicated devices, for example, one ormore field programmable gated array (FPGA) devices, to implement one ormore of the engines.

In an exemplary embodiment, the engines 152, 160, 164, 170, 172, 174,176, 177, 178, and 182 may be implemented with one or more dedicateddevices, for example, one or more field programmable gated array (FPGA)devices.

In an exemplary embodiment, the engines 152, 160, 164, 170, 172, 174,176, 177, 178, and 182 may be implemented with one or more softwaremodules comprising one or more software instructions stored on one ormore computer-readable mediums and executable by a processing unithaving one or more processors to cause the processing unit to performthe functions of the engines.

In an exemplary embodiment, the engines 152, 160, 164, 170, 172, 174,176, 177, 178, and 182 may be implemented with a combination of one ormore dedicated devices and/or one or more software modules stored on oneor more computer-readable mediums.

Embodiments described above may be implemented on the same computer,computer system, or a processor. Embodiments described above may also beimplemented on different computers, computer systems or processors.Further, embodiments described above may be implemented with one or moresoftware instructions or code residing on tangible computer-readablemedium or mediums of one or more computers or computer systems.

No element, act, or instruction used in the description of the exemplaryembodiments should be construed as critical or essential to theexemplary embodiments unless explicitly described as such.

The invention is described in detail with respect to exemplaryembodiments, and it will now be apparent from the foregoing to thoseskilled in the art that changes and modifications may be made withoutdeparting from the invention in its broader aspects, and the invention,therefore, as defined in the claims is intended to cover all suchchanges and modifications as fall within the true spirit of theinvention.

1. A magnetic resonance imaging apparatus, comprising: a magneticresonance imaging scanner to acquire diffusion-weighted imaging datafrom a subject; a reconstruction engine to reconstruct the acquireddiffusion-weighted imaging data into diffusion-weighted imagerepresentations of the subject; a diffusion tensor engine to construct adiffusion tensor map of an area of interest of the subject from thediffusion-weighted image representations, the diffusion tensor mapcomprising at least one voxel of tensor; an eigenvalue/eigenvectorordering engine to obtain and order pairs of eigenvectors andeigenvalues for at least one voxel in the constructed diffusion tensormap of the area of interest of the subject; a covariance matrixdetermining engine to construct a covariance matrix of a majoreigenvector for at least one voxel, with pairs of eigenvectors andeigenvalues obtained and ordered, in the constructed diffusion tensormap of the area of interest of the subject; a first normalized measuredetermining engine to compute a first normalized uncertainty measure ofat least one voxel, for the constructed covariance matrix of the majoreigenvector in the constructed diffusion tensor map of the area ofinterest of the subject; a second normalized measure determining engineto compute a second normalized uncertainty measure of at least onevoxel, for the constructed covariance matrix of the major eigenvector inthe constructed diffusion tensor map of the area of interest of thesubject; and a rendering engine to generate a human-viewable display ofan image representation of uncertainty from at least one of the firstnormalized uncertainty measure or the second normalized uncertaintymeasure.
 2. The magnetic resonance imaging apparatus according to claim1, wherein the first normalized uncertainty measure is an arealnormalized measure based on a ratio of a region on a unit sphereenclosed by a simple closed curve to an area of a hemisphere.
 3. Themagnetic resonance imaging apparatus according to claim 2, wherein theareal normalized measure is:${\Gamma = {\frac{2\; a}{\pi\; b\sqrt{1 + a^{2}}}\left\lbrack {{\left( {1 + b^{2}} \right){\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)}} - {K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}} \right\rbrack}},$where a and b are radii of an ellipse projected into a plane, and K andΠ are, respectively, the complete elliptic integrals of the first kindand the third kind.
 4. The magnetic resonance imaging apparatusaccording to claim 2, further comprising: a normalized areal mapconstructing engine to construct an areal map of the area of interest ofthe subject based on the areal normalized measure.
 5. The magneticresonance imaging apparatus according to claim 1, wherein the secondnormalized uncertainty measure is a normalized circumference based on aratio of a circumference of a simple closed curve on a unit sphere to acircumference of a great circle of the unit sphere.
 6. The magneticresonance imaging apparatus according to claim 5, wherein the normalizedcircumference is computed as:${\Lambda\left( {a,b} \right)} = {\frac{2}{\pi\; b\sqrt{1 + a^{2}}}\begin{pmatrix}{{\left( {1 + b^{2}} \right){\Pi\left( {\frac{a^{2} - b^{2}}{1 + a^{2}},\frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}}} \right)}} -} \\{K\left( \frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}} \right)}\end{pmatrix}}$ where a and b are radii of an ellipse of an ellipseprojected into a plane, and K and Π are, respectively, the completeelliptic integrals of the first kind and the third kind.
 7. The magneticresonance imaging apparatus according to claim 5, further comprising: anormalized circumference map constructing engine to construct anormalized circumference map of the area of interest of the subjectbased on the normalized circumference.
 8. The magnetic resonance imagingapparatus according to claim 1, further comprising: a cone ofuncertainty constructing engine to determine a cone of uncertaintysurrounding the major eigenvector for at least one voxel in theconstructed diffusion tensor map of the area of interest of the subjectbased on the determined first and second normalized uncertaintymeasures, the determined cone of uncertainty to be used by the renderingengine to generate as a human-viewable display of an imagerepresentation of uncertainty.
 9. A magnetic resonance imaging method,comprising: scanning a subject with a magnetic resonance imaging scannerto acquire diffusion-weighted imaging data; acquiring athree-dimensional diffusion tensor map of an area of interest of thesubject based on the acquired diffusion-weighted imaging data;processing the diffusion tensor map at a plurality of voxels in tensormap of the area of interest of the subject to obtain eigenvectors andeigenvalues; constructing covariance matrices of major eigenvectors forat least a portion of the plurality of voxels in the diffusion tensormap of the area of interest of the subject; computing first and secondnormalized uncertainty measures based on the constructed covariancematrices of the major eigenvectors for at least a portion of theplurality of voxels in the diffusion tensor map of the area of interestof the subject; and generating a human-viewable display of an imagerepresentation based on at least one of the first normalized uncertaintymeasures or the second normalized uncertainty measures for at least aportion of the plurality of voxels in the diffusion tensor map of thearea of interest of the subject.
 10. The method according to claim 9,wherein computing the first normalized uncertainty measures comprises:computing ratios of regions on unit spheres enclosed by simple closedcurves to areas of hemispheres.
 11. The method according to claim 10,further comprising: computing the ratios of regions on unit spheresenclosed by simple closed curves to areas of hemispheres as${\Gamma = {\frac{2\; a}{\pi\; b\sqrt{1 + a^{2}}}\left\lbrack {{\left( {1 + b^{2}} \right){\Pi\left( {{- b^{2}},\frac{a^{2} - b^{2}}{1 + a^{2}}} \right)}} - {K\left( \frac{a^{2} - b^{2}}{1 + a^{2}} \right)}} \right\rbrack}},$where a and b are radii of an ellipse projected into a plane, and K andΠ are, respectively, the complete elliptic integrals of the first kindand the third kind.
 12. The method according to claim 10, furthercomprising: constructing areal maps of the area of interest of thesubject based on the first normalized uncertainty measures.
 13. Themethod according to claim 9, wherein computing the second normalizeduncertainty measure comprises: computing ratios of circumferences ofsimple closed curves on unit spheres to circumferences of great circlesof the unit spheres.
 14. The method according to claim 13, furthercomprising: computing the ratios of circumferences of simple closedcurves on unit spheres to circumferences of great circles of the unitspheres as:${\Lambda\left( {a,b} \right)} = {\frac{2}{\pi\; b\sqrt{1 + a^{2}}}\begin{pmatrix}{{\left( {1 + b^{2}} \right){\Pi\left( {\frac{a^{2} - b^{2}}{1 + a^{2}},\frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}}} \right)}} -} \\{K\left( \frac{b^{2} - a^{2}}{\left( {1 + a^{2}} \right)b^{2}} \right)}\end{pmatrix}}$ where a and b are radii of an ellipse of an ellipseprojected into a plane, and K and Π are, respectively, the completeelliptic integrals of the first kind and the third kind.
 15. The methodaccording to claim 13, further comprising: constructing normalizedcircumference maps of the area of interest of the subject based on thesecond normalized uncertainty measures.
 16. The method according toclaim 9, further comprising: determining cones of uncertaintysurrounding the major eigenvectors for the plurality of voxels based onthe first and second normalized uncertainty measures, the determinedcones of uncertainty to be used for generating a human-viewable displayof an image representation of uncertainty.
 17. A magnetic resonanceimaging apparatus, comprising: means for acquiring diffusion-weightedimaging data from a subject; means for reconstructing the acquireddiffusion-weighted imaging data into diffusion-weighted imagerepresentations of the subject; means for constructing a diffusiontensor map of an area of interest of the subject from the reconstructeddiffusion-weighted image representations, said diffusion tensor mapcomprising at least one voxel of tensor; means for obtaining andordering pairs of eigenvectors and eigenvalues for at least one voxel inthe constructed diffusion tensor map of the area of interest of thesubject; means for constructing a covariance matrix of a majoreigenvector for at least one voxel, with pairs of eigenvectors andeigenvalues obtained and ordered, in the constructed tensor map of thearea of interest of the subject; means for computing a first normalizeduncertainty measure of at least one voxel, for which the covariancematrix of the major eigenvector has been constructed, in the constructedtensor map of the area of interest of the subject; means for computing asecond normalized uncertainty measure of at least one voxel, for whichthe covariance matrix of the major eigenvector has been constructed, inthe constructed tensor map of the area of interest of the subject; andmeans for generating a human-viewable display of an image representationof uncertainty based on at least one of the first normalized uncertaintymeasure or the second normalized uncertainty measure.
 18. One or morecomputer-readable mediums comprising software, which software, whenexecuted by a computer system, causes the computer to perform operationsto analyze diffusion-weighted imaging representations reconstructed fromacquired diffusion-weighted imaging data of a subject, thecomputer-readable medium comprising: one or more instructions executableby the computer to construct a diffusion tensor map of an area ofinterest of the subject from the reconstructed diffusion-weighted imagerepresentations of the subject, said diffusion tensor map comprising atleast one voxel of tensor; one or more instructions executable by thecomputer to obtain and order pairs of eigenvectors and eigenvalues forat least one voxel in the constructed diffusion tensor map of the areaof interest of the subject; one or more instructions executable by thecomputer to construct a covariance matrix of a major eigenvector for atleast one voxel, with pairs of eigenvectors and eigenvalues obtained andordered, in the constructed diffusion tensor map of the area of interestof the subject; one or more instructions executable by the computer tocompute a first normalized uncertainty measure of at least one voxel,for which the covariance matrix of the major eigenvector has beenconstructed, in the constructed diffusion tensor map of the area ofinterest of the subject; one or more instructions executable by thecomputer to compute a second normalized uncertainty measure of at leastone voxel, for which the covariance matrix of the major eigenvector hasbeen constructed, in the constructed diffusion tensor map of the area ofinterest of the subject; and one or more instructions executable by thecomputer for generating a human-viewable display of an imagerepresentation of uncertainty based on at least one of the firstnormalized uncertainty measure.